/*
 *	Copyright (C) 2008 CRIMERE
 *	Copyright (C) 2008 Jean-Marc Mercier
 *	
 *	This file is part of OTS (optimally transported schemes), an open-source library
 *	dedicated to scientific computing. http://code.google.com/p/optimally-transported-schemes/
 *
 *	CRIMERE makes no representations about the suitability of this
 *	software for any purpose. It is provided "as is" without express or
 *	implied warranty.
 *
 *  This program is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation, either version 3 of the License, or
 *  (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *   You should have received a copy of the GNU General Public License
 *   along with this program.  If not, see <http://www.gnu.org/licenses/>.
 *
 */
#include <set>


#include <src/math/Operators/u_2_Delta.h>
#include <src/algorithm/Evolution_Algorithm/Crank_Nicolson.h>
#include <src/math/CauchyProblems/CauchyProblem_traits.h>


//! Cauchy Problem for the non linear heat equation \f$\partial_t u(t,x) = u^2 \Delta u(t,x)\f$.
/*!
	This class initialize the Cauchy problem for the one dimensional non linear heat equation
	\f[
		\partial_t u(t,x) = u^2 \Delta u(t,x)
	\f];
	with \f$x \in [0,1]\f$, \f$t\ge 0\f$.

	It illustrate the behavior of the finite difference scheme
	\f[
		\delta_t u^n = \left(u^{n+1/2}\right)^2\delta_{x\bar{x}} u^{n+1/2} = 0
	\f];
	The algorithm used is Crank Nicolson.
	The boundary conditions are given by <data> : either Dirichlet or Periodic ones.
	The initial data corresponds to a normal distribution centered in 0.5.
*/
#define NUMBER_OF_PARTICLES 20

#if !defined(_Cauchy_u_2_Delta_)
#define _Cauchy_u_2_Delta_


template<class data>
struct CrankNicolson_u_2_Delta : public Crank_Nicolson<u_2_Delta_t<data>>
{
//	OS_STATIC_CHECK(data::Dim == 1); // not sure of the multidimensional behavior. To test before.
	virtual OS_double CFL() {
		return positive_define_cfl<time_operator_type>()(iterate, Op_iterate);
	};
	virtual bool criteria() {
		OS_double breakpoint1 = jmmath::norm22(iterate - S_k1);
		OS_double breakpoint2 = jmmath::norm22(iterate);
		return breakpoint1/(1.+breakpoint2) < 1e-4;
	};
};

template <class data, std::size_t Dim, class Interface = Graphical_Interface>
class Cauchy_u_2_Delta : public CauchyProblem_traits< CrankNicolson_u_2_Delta < data > , Graphical_Interface	>
{
	OS_STATIC_CHECK(false); // Test before validating.
};

template <class data>
		class Cauchy_u_2_Delta<data,2,Graphical_Interface> : public CauchyProblem_traits< CrankNicolson_u_2_Delta < data > , Graphical_Interface	>
		{
		public :
			typedef typename public CauchyProblem_traits< Crank_Nicolson < u_2_Delta_t<data> > , Graphical_Interface	> Base;
			Cauchy_u_2_Delta(){
				initial_conditions(NUMBER_OF_PARTICLES);
			};
			virtual smart_ptr<interface_type> Create_Instance() {
				return smart_ptr<interface_type>(new Cauchy_u_2_Delta);
			};

			virtual void set_dim(OS_size dim) {
				OS_DYNAMIC_CHECK(dim == 2,"Cauchy_u_2_Delta: only dimension 2");
			};
			OS_size get_dim() {return 1;};
			virtual void initial_conditions(OS_size nb_part)
			{
				std::vector<std::size_t> size_(data::Dim);
				std::size_t inf = std::pow(nb_part,1./((double) data::Dim));
				std::fill(size_.begin(),size_.end(),inf+1);
				std::vector<double> average(data::Dim);
				std::fill(average.begin(),average.end(),.5);
				set_initial_conditions(Gaussian_Functor<result_type>(average)(size_));
			};
			virtual typename interface_type::p_array_2 get_representation() {
				return representation_3<data,2,Finite_Diff_Operator>()(get_solution());
			};
		};

template <class data>
		class Cauchy_u_2_Delta<data,1,Graphical_Interface> : public CauchyProblem_traits< CrankNicolson_u_2_Delta < data > , Graphical_Interface	>
		{
		public :
			typedef typename public CauchyProblem_traits< Crank_Nicolson < u_2_Delta_t<data> > , Graphical_Interface	> Base;
			Cauchy_u_2_Delta(){
				initial_conditions(NUMBER_OF_PARTICLES);
			};
			virtual smart_ptr<interface_type> Create_Instance() {
				return smart_ptr<interface_type>(new Cauchy_u_2_Delta);
			};

			virtual void set_dim(OS_size dim) {
				OS_DYNAMIC_CHECK(dim == 1,"Cauchy_u_2_Delta: only dimension 1");
			};
			OS_size get_dim() {return 1;};
			void initial_conditions(OS_size nb_part)
			{
				set_initial_conditions(Gaussian_Functor<data>()(nb_part));
			};
		};

		
//! Cauchy Problem for the one dimensional heat equation \f$\partial_t u = \partial_{xx} u\f$, particle method.
/* 
    This class initialize the Cauchy problem for the one dimensional heat equation
    \f[
        \partial_t u(t,x) = \partial_{xx} u(t,x)
    \f]
    with \f$x \in \RR \f$, \f$t\ge 0\f$, \f$u(0,x) = u_0(x)>0\f$, \f$\int_{\RR} u_0 = 1\f$. We solve this equation by the particle method,
    considering its Lagrangian form
    \f[
        \partial_t S = -\partial_y \frac{1}{\partial_y S}
    \f]
    with \f$y \in ]0,1[ \f$, \f$t\ge 0\f$, with the boundary conditions \f$(S(t,0)=-\infty;S(t,1)=+\infty)\f$. \f$S(t,y)\f$ is the increasing rearrangement of \f$u(t,x)\f$ : \f$u(S(t,y))=\frac{1}{\partial_y S(t,y)}\f$. This equation may be solved directly by a
    direct method (see the Cauchy Problem ) . We propose here a different method to solve this last equation : let \f$u(t,y)= u(S(t,y))=\frac{1}{\partial_y S(t,y)}\f$. We compute that
    \f$u(y)\f$ verifies the following equation :
    \f[
        \partial_t u(y) = u^2(y)\partial_{yy} u(t,y)
    \f]
    with \f$y \in [0,1] \f$, \f$t\ge 0\f$, and Dirichlet boundary conditions \f$u(0)=u(1)=0\f$. This equation is interesting from a numerical point of view because the operator \f$u^2\partial_{yy} u\f$
    is positive defined : i.e. for all \f$\varphi \in \mathcal{D}^{0,1}([0,1])\f$
    \f[
        \varphi^2\partial_{yy} \varphi \cdot \varphi = \int_{[0,1]} \varphi^3\partial_{yy}\varphi  \le 0
    \f]
    We use the Crank Nicolson scheme to compute this equation :
    \f[ \label{Delta_u_2}
        \delta_t u^n = (u^{n+1/2})^2 \delta_{y\bar{y}} u^{n+1/2}
    \f]
    This scheme verifies the discrete entropy conservation
    \f[
        \delta_t U(u^n) = \delta_t \sum_i (u_i^n)^2 = \sum_i (u_i^{n+1/2})^3 \delta_{x\bar{x}} u_i^{n+1/2} \le 0
    \f]
    Its natural CFL is the one given by positive defined operators : denoting \f$\mathcal{O}(S^{n})= (u^{n})^2 \delta_{y\bar{y}} u^{n}\f$
we define the n-th time step to be (see positive define cfl)
    \f[
        \tau^n =  -2 \frac { \mathcal{O}(S^{n}) \cdot S^{n} } {|\mathcal{O}(S^{n})|_{\ell^2}^2}
    \f]
    The solution is represented as \f$\{S_i^n,u_i^n\}_{i=0,..,N;n\ge 0}\f$, where the particle location
    location \f$\{S_i^n\}_i\f$ may be found either by inverting the system of equations
    \f[
        u_i^n = \frac{1}{\delta_{y} S_i^n}
    \f]
    or solving the equation
    \f[
        \delta_t S_i^n = -\delta_{\bar{y}} \frac{1}{\delta_{y} S_i^n}
    \f]
*/
template <class data>
class Cauchy_Delta_p_change_1 : public Cauchy_u_2_Delta< data, data::Dim,Graphical_Interface>
		{
		OS_STATIC_CHECK(data::Dim ==1);
		public :
			Cauchy_Delta_p_change_1(){
				initial_conditions(NUMBER_OF_PARTICLES);
			};

			virtual smart_ptr<interface_type> Create_Instance() {
				return smart_ptr<interface_type>(new Cauchy_Delta_p_change_1);
			};

			virtual void set_dim(OS_size dim) {
				OS_DYNAMIC_CHECK(dim == 1,"Not Yet Implemented");
			};
			OS_size get_dim() {return 1;};
			void initial_conditions(OS_size nb_part)
			{
				get_graph_.set_state_variable(smart_ptr<data> (new data(nb_part)) ); // for the graph
				set_initial_conditions(particle_density<data>() ( InverseCumulativeFunctor<result_type>()(nb_part)) );
			};
			virtual typename interface_type::p_array_2 get_representation() {
			// we cheat because we are lazy : we compute explicitely the particle location through the inverse normal function.
				const data& solution_ = get_solution();
				data p_ (InverseCumulativeFunctor<data>(0.,get_current_time() + std::sqrt( 0.1 )) (solution_.size()) );

				p_array_2 temp(solution_.size(),2);
				for (OS_size i = 0; i < temp.size(); ++i) {
					temp[i][0] = p_[i];
					temp[i][1] = solution_[i];
				};
				return temp;
			};
			virtual std::string get_name() {return "Particle 1D heat Equation with variable change."
													+time_operator_type::get_name()
													+time_operator_type::Boundary_type::get_name()
													+ result_type::get_name()
													+". "
													+ time_operator_type::Family_Operator::get_name();
			};
			virtual const graph_type& get_graph() const
			{
				return get_graph_.get_graph();
			};

			u_2_Delta_t<data> get_graph_;

		};

template<class data>
struct Cauchy_Strange_interface : public Crank_Nicolson<Delta_t<data,data::Dim>>
{
	typedef Crank_Nicolson<Delta_t<data,data::Dim>> Base;
	typedef typename data::base_value_type  base_value_type; 
	typedef Neumann_Vector<Neumann_Vector<base_value_type>> neuman_type;
	typedef OS_Vector<neuman_type> S_type;
	Cauchy_Strange_interface() {
		OS_STATIC_CHECK(data::Dim == 2); // "Cauchy_Delta_p_change_interface : not sure of the multidimensional behavior. Test and validate before.");
	};
	virtual OS_double CFL() {
		S_12 = S_;
		S_iterate_ = S_;
		S_k12 = *get_state_variable();
		breakpoint2 = jmmath::norm22(*get_state_variable());
		return positive_define_cfl<time_operator_type>()(S_k12, Op_(t_12,S_k12)^pow(S_k12,2));
	};
	virtual void vstar() {
		S_k12   = (*get_state_variable()+iterate)/2.;
		S_12    = (S_iterate_ + S_)/2.;
	};
	virtual bool criteria() {
		OS_double breakpoint1 = jmmath::norm22(iterate - S_k1);
		if (breakpoint1/(1.+breakpoint2) < 1e-4) // this end the Crank_Nicolson loop
		{
			S_ = S_iterate_; // update positions
			return true;
		};
		return false;
	};
	virtual void update(){
		Op_.set_state_variable(get_state_variable());
		grad_c_.set_state_variable(get_state_variable());
		if (S_.capacity(data::Dim) != get_state_variable()->capacity()){
			typedef  S_type::indice_list_type indice_list_type; 
			indice_list_type local(S_type::Dim);
			indice_list_type indx(get_state_variable()->size_list());
			local[S_type::Dim-1] = data::Dim;
			std::copy(indx.begin(),indx.end(),local.begin());
			S_.resize(local);	
		};
	};
	virtual void gstar() {
		S_iterate_ = S_ - (delta_t_*data::Dim)*(grad_c_( S_k12));
		Op_iterate = data::Dim*( S_k12^S_k12^(Op_(t_12,S_k12)) );
	};
	virtual void adjust_cfl(const time_type & from, argument_type S_k){
		S_12 = S_;
		S_iterate_ = S_;
		Base::adjust_cfl(from,S_k);
	}
	S_type S_, S_iterate_, S_12;
	Grad_c_<argument_base_type,S_type> grad_c_;
	OS_double breakpoint2;
};
/*! \brief Cauchy Problem for the two dimensional equation \f$\partial_y S = - D \nabla \frac{1}{\nabla \cdot S}\f$.

Note : the following analysis is multidimensional. The method is implemented only for Dim==2 at present time.
We study here the following equation :
    \f[ \label{whynot}
        \partial_y S = - D \nabla \frac{1}{\nabla \cdot S}
    \f]
where \f$S \equiv S(t,y)\f$, \f$y=\{y_i\}_{i=1,..,D} \in [0,1]^D\f$, and
initial conditions \f$S(0,y) = S_0(y)\f$. This equation has some
interesting properties : let us suppose that \f$S_0 \in
\mathcal{H}^1([0,1]^D)\f$ and \f$\nabla \cdot S_0
> 0\f$, then
    \f[
        \frac{1}{2}\frac{d}{dt}\int_{[0,1]^D} \|S\|^2 = - D \int_{[0,1]^D} S \cdot \nabla \frac{1}{\nabla \cdot
        S} = D \int_{[0,1]^D}\nabla \cdot S \frac{1}{\nabla \cdot S} = D.
    \f]
Let us define the map \f$x=S(y)\f$, and denote \f$u(t,x) = |\nabla S(y)|\f$ the determinant of the Jacobian of \f$S\f$. Equivalently, the previous relation gives us
    \f[
     \frac{d}{dt} \int_{\RR^D}\|x\|^2u(t,x)dx = \frac{d}{dt} \int_{\RR^D}\|S\|^2dy = 2D
    \f]
that is characteristic of the heat equation. We don't know whether this map determines a transport map or not. 
Any suggestion are welcomed !


For the numerical simulation, we will be considering the following system, equivalent to (\ref{whynot})
 \f[
\label{system}
    \begin{split}
    \partial_t S = - \nabla v\\
\partial_t v = D v^2 \Delta v
    \end{split}
\f]

Submitted to initial conditions \f$\{S_0,u_0\}\f$. Let us determine initial condition corresponding to an
optimal transport map : let \f$\sigma >0\f$ a real and consider the gaussian function
    \f[
        \mathcal{N}_{\sigma}(x) = \frac{1}{\sqrt{2 \pi}\sigma} \exp\left(-\frac{x^2}{\sigma}\right)
    \f]
    for \f$x\in \RR\f$. Let us call its integral \f$s^{-1}(x)=\int_{-\infty}^x \mathcal{N}_{\sigma}(u)
    du\f$. We recall that \f$s^{-1}(x) \in [0,1]\f$ transports optimally
    \f$\mathcal{N}_{\sigma}\f$ into the uniform distribution of the
    segment \f$[0,1]\f$. Let us call now \f$s(y)\f$ its inverse mapping, that
    transports the uniform distribution of \f$[0,1]\f$ into
    \f$\mathcal{N}_{\sigma}\f$. Since \f$s\f$ is an optimal transport map,
    there exists a convex function \f$\phi(y)\f$ such that \f$s= \partial_y
    \phi\f$. We recall also that, derivating the trivial relation
    \f$s^{-1}(s(y))=y\f$, we deduce
    \f[
        \partial_y s = \frac{1}{\mathcal{N}_{\sigma}(s)}
    \f]


    Now let us define the \f$D\f$ dimensional function of the hypercube
    \f$y=\{y_i\}_{i=1,..,D} \in [0,1]^D\f$
    \f[
        \Phi(y) := \sum_{i=1,..,D} \phi(y_i)
    \f]
    \f$\Phi\f$ is clearly a convex function. Thus it represents an
    optimal transport map that we call \f$S(y)\f$, with
    \f[
        S(y) = \nabla \Phi = \{ s(y_i) \}_{i=1,..,D}
    \f]
    We compute also that
    \f[
        \nabla \cdot S(y) = \sum_{i=1,..,D} \frac{1}{\mathcal{N}_{\sigma}(y_i)}
    \f]
    Thus we choose for the initial conditions of our system
    (\ref{system}) the following ones :
    \f[
        S_0(y) = \sum_{i=1,..,D} s(y_i); u_0(y) = \frac{1}{\sum_{i=1,..,D} \frac{1}{\mathcal{N}_{\sigma}(y_i)}}
    \f]
*/
template <class data, class A = Cauchy_Strange_interface<data> > //default is heat equation.
		class Cauchy_Strange : public CauchyProblem_traits< A , Graphical_Interface	>
		{
		public :
			typedef Cauchy_Strange type;
			typedef CauchyProblem_traits Base;

			type(){
				initial_conditions(NUMBER_OF_PARTICLES);
			};
			virtual smart_ptr<interface_type> Create_Instance() {
				return smart_ptr<interface_type>(new type);
			};
			virtual void initial_conditions(OS_size nb_part)
			{
				std::size_t inf = std::pow(nb_part,1./((double) data::Dim));
				data initial_conditions_(inf+1,inf+1);

				typedef Unbounded_Vector<data::base_value_type> unbounded_array_1;

				InverseCumulativeFunctor<unbounded_array_1> i_g;
				unbounded_array_1 position(i_g(inf+1)); // s

				OS_normal_distribution<> gaussian_function(i_g.function_.average(), i_g.function_.sigma());
				data::value_type local(gaussian_function.operator()<data::value_type,data::value_type>(position));
				for (data::indice_type i=0; i < inf+1 ; ++i)
					for (data::indice_type j=0; j < inf+1 ; ++j)
						initial_conditions_[i][j] = 1./ (1./local[i] + 1./local[j]);

				set_initial_conditions(initial_conditions_);

				typedef algorithm_type::S_type S_type;
				S_type& S = get_algorithm().S_;
				for (data::indice_type i=0; i < inf+1 ; ++i) {
						for (data::indice_type j=0; j < inf+1 ; ++j) {
							S[0][i][j] = position[i];
							S[1][i][j] = position[j];
						};
					};
						
			};
			virtual std::string get_name() {return "2D partial_t S = - D nabla frac{1}{ nabla cdot S}";};

			virtual void set_dim(OS_size dim) {
				OS_DYNAMIC_CHECK(dim == 2,"Cauchy_Strange: only dimension 2");
			};
			OS_size get_dim() {return 2;};
			virtual typename interface_type::p_array_2 get_representation() {
				typedef algorithm_type::S_type S_type;
				S_type& S = get_algorithm().S_;
				const data& v = get_solution();

				OS_size size_1 = v.size(1);
				OS_size size_2 = v.size(2);
				OS_double alpha = .1;
				OS_double coef1 = std::sin(alpha);
				OS_double coef2 = std::cos(alpha);
				p_array_2 temp(size_1*size_2,2);
				for (OS_size i = 0; i < size_1; ++i) {
					for (OS_size j = 0; j < size_2; ++j) {
						double breakpoint = v[i][j];
						temp[i+j*size_1][0] = S[0][i][j]*coef2 - S[1][i][j]*coef1;
						temp[i+j*size_1][1] = S[1][i][j]*coef1 + breakpoint;
//						temp[i+j*size_1][0] = S[0][i][j];
//						temp[i+j*size_1][1] = S[1][i][j];
					};
				};
				return temp;
			};
		};

template<class data>
struct Cauchy_Delta_p_change_interface : public Crank_Nicolson<Delta_t<data,data::Dim>>
{
	typedef Crank_Nicolson<Delta_t<data,data::Dim>> Base;
	typedef typename data::base_value_type  base_value_type; 
	typedef typename data::value_type  data_value_type; 
	typedef Neumann_Vector<Neumann_Vector<base_value_type>> neuman_type;
	typedef OS_Vector<neuman_type> S_type;
	Cauchy_Delta_p_change_interface() {
		OS_STATIC_CHECK(data::Dim == 3); // "Cauchy_Delta_p_change_interface : not sure of the multidimensional behavior. Test and validate before.");
	};
	virtual OS_double CFL() {
		S_12 = S_;
		S_iterate_ = S_;
		S_k12 = *get_state_variable();
		breakpoint2 = jmmath::norm22(*get_state_variable());
		return positive_define_cfl<time_operator_type>()( S_k12, Op_iterate^S_k12^S_k12 );
	};
	virtual void vstar() {
		OS_Array_Copy(S_k12,(*get_state_variable()+iterate)/2.);
		OS_Array_Copy(S_12, (S_iterate_ + S_)/2.);
	};
	virtual bool criteria() {
		OS_double breakpoint1 = jmmath::norm22(iterate - S_k1);
		if (breakpoint1/(1.+breakpoint2) < 1e-4) // this end the Crank_Nicolson loop
		{
			S_ = S_iterate_; // update positions
			return true;
		};
		return false;
	};
	virtual void update(){
		Op_.set_state_variable(get_state_variable());
		if (S_.capacity() != get_state_variable()->capacity()){
			typedef  S_type::indice_list_type indice_list_type; 
			indice_list_type local(get_state_variable()->size_list());
			S_.resize(local);
			OS_Array_Copy(S_,*get_state_variable());
			local.pop_back();
			u_1_u_2_.resize(local);
			S_iterate_.resize(local);
			Op_iterate.resize(local);
		};
	};
	virtual void gstar() {
		S_k12[0]^S_k12[1];
		OS_Array_Copy(u_1_u_2_,S_k12[0]^S_k12[1]);
		OS_Array_Copy(S_iterate_ ,S_ - (delta_t_)*(grad_c_(S_k12[0]+S_k12[1])));
		OS_Array_Copy(Op_iterate,S_k12^S_k12^(Op_(t_12,S_k12)));
	};
	virtual void adjust_cfl(const time_type & from, argument_type S_k){
		S_12 = S_;
		S_iterate_ = S_;
		Base::adjust_cfl(from,S_k);
	}
	S_type S_, S_iterate_, S_12;
	data_value_type u_1_u_2_;
	OS_double breakpoint2;
	Grad_c_<data_value_type,S_type,data_value_type::Dim,smart_ptr<S_type> > grad_c_;
};
template <class data, class A = Cauchy_Delta_p_change_interface<data> > //default is heat equation.
		class Cauchy_Delta_p_change_2 : public CauchyProblem_traits< A , Graphical_Interface	>
		{
		public :
			typedef CauchyProblem_traits Base;

			Cauchy_Delta_p_change_2(){
				initial_conditions(NUMBER_OF_PARTICLES);
			};
			virtual smart_ptr<interface_type> Create_Instance() {
				return smart_ptr<interface_type>(new Cauchy_Delta_p_change_2);
			};
			virtual void initial_conditions(OS_size nb_part)
			{

				std::size_t inf = std::pow(nb_part,1./((double) data::Dim-1.));

				data initial_conditions_(data::Dim-1,inf+1,inf+1);

				typedef Unbounded_Vector<data::base_value_type> unbounded;
				typedef Dirichlet_Vector<data::base_value_type> dirichlet;
				InverseCumulativeFunctor<unbounded> i_g;
				unbounded position(i_g(inf+1)); // s

				Gaussian_Function<> gaussian_function(i_g.function_.average(), i_g.function_.sigma());
				dirichlet local(gaussian_function.operator()<dirichlet,dirichlet>(position));

				for (data::indice_type i=0; i < inf+1 ; ++i)
					for (data::indice_type j=0; j < inf+1 ; ++j) {
						initial_conditions_[0][i][j] = local[i];
						initial_conditions_[1][i][j] = local[j];
					};
				set_initial_conditions(initial_conditions_);
						
				typedef algorithm_type::S_type S_type;
				S_type& S = get_algorithm().S_;
				for (data::indice_type i=0; i <inf+1 ; ++i) {
						for (data::indice_type j=0; j < inf+1 ; ++j) {
							S[0][i][j] = position[i];
							S[1][i][j] = position[j];
						};
					};

			};
			virtual std::string get_name() {return "2D partial_t u = Delta u. Particle method with variable change.";};

			virtual void set_dim(OS_size dim) {
				OS_DYNAMIC_CHECK(dim == 2,"Cauchy_u_2_Delta: only dimension 2");
			};
			OS_size get_dim() {return 2;};
			virtual typename interface_type::p_array_2 get_representation() {
				typedef algorithm_type::S_type S_type;
				S_type& S = get_algorithm().S_;
				const data& v = get_solution();

				OS_size size_1 = v.size(1);
				OS_size size_2 = v.size(2);
				OS_double alpha = .1;
				OS_double coef1 = std::sin(alpha);
				OS_double coef2 = std::cos(alpha);
				p_array_2 temp(size_1*size_2,2);
				for (OS_size i = 0; i < size_1; ++i) {
					for (OS_size j = 0; j < size_2; ++j) {
						double breakpoint = v[0][i][j]*v[1][i][j];
						temp[i+j*size_1][0] = S[0][i][j]*coef2 - S[1][i][j]*coef1;
						temp[i+j*size_1][1] = S[1][i][j]*coef1 + breakpoint;
//						temp[i+j*size_1][0] = S[0][i][j];
//						temp[i+j*size_1][1] = S[1][i][j];
					};
				};
				return temp;
			};
		};


#undef NUMBER_OF_PARTICLES

#endif